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ABSTRACT 

Due to the gravitational influence of density fluctuations driven by magneto- 
rotational instability in the gas disk, planetesimals and protoplanets undergo diffusive 
radial migration as well as changes in other orbital properties. The magnitude of the 
effect on particle orbits can have important consequences for planet formation scenarios. 
We use the local-shearing-box approximation to simulate an ideal, isothermal, magne- 
tized gas disk with vertical density stratiflcation and simultaneously evolve numerous 
massless particles moving under the gravitational field of the gas and the host star. We 
measure the evolution of the particle orbital properties, including mean radius, eccen- 
tricity, inclination, and velocity dispersion, and its dependence on the disk properties 
and the particle initial conditions. Although the results converge with resolution for 
fixed box dimensions, we find the response of the particles to the gravity of the turbulent 
gas correlates with the horizontal box size, up to 16 disk scale heights. This correlation 
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indicates that caution should be exercised when interpreting local-shearing-box models 
involving gravitational physics of magneto-rotational turbulence. Based on heuristic 
arguments, nevertheless, the criterion L^/R ~ 0(1), where Lh is the horizontal box size 
and R is the distance to the host star, is proposed to possibly circumvent this conun- 
drum. If this criterion holds, we can still conclude that magneto-rotational turbulence 
seems likely to be ineffective at driving either diffusive migration or collisional erosion 
under most circumstances. 

Subject headings: Instabilities — Magnetohydrodynamics (MHD) — Planets and satel- 
lites: formation — Planet-disk interactions — Protoplanetary disks — Turbulence 



1. INTRODUCTION 



The gas disk surrounding a young stellar object is usually argued to be a turbulent accretio n 
disk when the star is still in the process of accumulating materials (see, e.g., iFrank et al.l |2002| ) . 
While the molecular viscosity of the gas is negligibly small, shear stress resulting from turbulence 
can effectively drive angular momentum transport in the disk. At the same time, the dynamics 
of protoplanetary objects embedded in the disk is significantly affected by their interactions with 
the turbulent gas, at a stag e when these objects are still amassing solids to grow in size (e.g., 
Papaloizou fc Terquemll2006l . and references therein). Understanding particle dynamics in a turbu- 
lent gas disk, therefore, will provide insight into the paths along which new planets may eventually 
form. 

One of the most promising mec hanisms to drive the t urbulence is through the magneto- 
rotational instability (MRI; see, e.g.. iBalbus &: HawlevI Il998l . and references therein). A weakly 



magnetized, differentially rotating gas disk is unstable to linear perturbations, and after nonlinear 
effects set in, the gas presumably reaches a statistically steady, sustained, turbulent state. The 
underlying processes that deter r nine the properties of this saturat e d state are still under active re- 
searc h (e.g., Ijamroz et ahlboosl : Latter et ahlbood : Ivishniadbood : IPessahlboiol : loishi fc Mac Low 



20 111 ). Nevertheless, it is known that the stren gth of this ma,gneto - rotational turbulence depends on 



the n et magnetic flux threading the disk (e.g.. iHawley et al.lll995l : IJohansen et al.ll2006l : lYang et al 



20091 . hereafter YMM09), and thus can be controlled to the desired level and treated as a parameter. 



The gravitational torques exerted on an embedded solid object due to density fluctuations 
in magneto-rotational t urbulence have been shown to be stochastic, rendering the rnotion of the 



object a randoi n walk (jLaughlin et al 



2007; YMM09; 



disk (Johnson et al 



20041 : iNelson fc Papaloizoij l2004l : iNelsonI 1200,'tI : lOishi et al 



Nelson &: Gressell l2010l') oii top of the ordered migrati on expected in a laminar 



20061 : lAdams k BlochI l2009l : iKetchum et al.l I2OIII ). No te that the ordered 



migr a tion rnay be either inw ard or outward in different regions of the disk (jPaardekooper et al 



2010l . I2OIII : iLyra et al.l I2OIOI ). In addition to inducing this radial diffusive migration, the same 



gravitational force drives the diffusion of orbital eccentricity of the solid objects, increasing their 
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velocity dispersion (|Ogihara et alJ l200?l : llda et alJ I2OO8I : YMM09; iNelson Gresseil l20im . The 
effect on tlie orbital inclination, however, has not been investigated yet, as that requires the inclusion 
of the vertical component of the gravity from the host star. Overall, the gravitational influence of 
the turbulent gas on particle orbital properties may play an important role in the course of planet 
formation. 

To assess the significance of this process in planet formation scenarios, numerical simulations 
simultaneously evolving magneto-rotational turbulence and particle dynamics can be conducted and 
the resulting particle orbital evolution can be statistically measured. However, a consistent picture 
of the actual magnitude of this effect has been elusive. Using the local-shearing-box simulations, 
we reported in YMM09 that the response of test particles to the gravity of the turbulent gas is 
systematically lower than what was previously reported in global disk models. We also found in 
the same study t hat the effect signi f icantl y depends on the horizontal box size of the numerical 
models. Recently, iNelson &: Gressell (j2010|) have found similar behavior in their local models and 
argued that large box size might be necessary to see convergence in the stochastic torque and 
forcing generated by the turbulence. 

Extending the local unstratified disk models of YMM09, we now consider disks with vertical 
density stratification by including the linearized vertical gravity from the central star. We examine 
whether a stratified disk model gives results consistent with a comparable unstratified model and 
whether convergence in particle orbital evolution, if any, can be achieved with local-shearing-box 
simulations. In Section [2l we describe in detail our numerical models. We report the properties 
of our modeled magneto-rotational turbulence in Section [3] and the statistical evolution of particle 
orbital properties in Section HI We analyze the gravitational field generated by the turbulence in 
Section [5] and discuss the convergence issues of our local models in Section [6l We conclude our 
discussion in Section [71 



2. NUMERICAL MODELING 

Directly extending the work we reported in YMM09, we continue to use the Pencil Cod^ll to 
model particles moving in magneto-rotational turbulence. We describe in detail the equations and 
the relevant numerical techniques we implemented in the code. 



2.1. Magnetohydrodynamics 



We use the local shea ring box approximation (e.g., Goldreich fc Lvnden-Bellll965 : Brandenburg et al 



I995I : iHawlev et al.lll995l ). A local shearing box is a small Cartesian box at a large distance R from 



^The Pencil Code is publicly available at http://code.google.eom/p/pencil-code/. 
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the host star such that the center of the box revolves at the Keplerian angular speed ^1k- The 
box is always oriented with the x-axis directed radially and the y-axis azimuthally. In contrast to 
YMM09, we include the linearized vertical gravity from the host star so that the disk is vertically 
stratified. We again impose a vertical, external magnetic field -Bgxt = B^xt^ to maintain a non-zero 
magnetic flux. The MHD equations then become 

3 

dtp - -^Koodyp + V ■ {pu) = fo, (1) 

3 1/1 

dtu Q.KxdyU + u ■ Vu = — Vp + ( 2il,KUyX VlKUxy — ^\zz 

2 p \ 2 

+ ^Jx (B + Bext)+/y, (2) 
3 3 

dtA - -VLKxdyA = -QxAyX + ux {B + Bext) + fn, (3) 

where p is the gas density, u is the gas velocity relative to the background shear flow, p is the 
gas pressure, J = V x B/po is the electric current density, B = V x A, and po is the vacuum 
permeability. The terms fo, fv^ are numerical dissipation terms, including both hyper- 

diffusion and shock diffusion, that are needed to stabilize the scheme. The reader is referred to 
YMM09 for their description. 

To close the system of equations ([I])-©, we assume an isothermal "equation of state" , p = c^p, 
where Cs is the isothermal speed of sound. This is an ideal limit of the energy equation, in which heat 
diffusion is effectively instantaneous. Gas buoyancy should thus be absent a nd the vertical mixing i n 



the resulting magneto-rotational turbu lence might be artificially enhanced (jBai fc: Goodmanll2009l ) 



Nevertheless, lOishi fc: Mac Lowl (j2009l ) showed that vertical motions are not completely suppressed 
even in the case of adiabatic gas, where no heat diffusion occurs, and the differences in the mean 
kinetic and the magnetic energies of a vertically stratified gas disk between the isothermal and the 
adiabatic limits may only amount to about a factor of two. 

We set up the gas density so that the gas is in vertical hydrostatic equilibrium initially: 

Poiz) = pm exp (^-Jp^ ' (4) 

where pm is the mid-plane gas density and H = \f2csl^K is the vertical disk scale height H We 
fix the vertical dimension of the computational domain at \z\ < 2H. The initial magnetic vector 
potential A is set to zero, while the external magnetic field is fixed at such a level that the corre- 
sponding plasma /3ext(-z) = '^Poc^Po/B^xt is /5ext(0) = 6.2 x 10^ in the mid-plane. Gaussian noise of 
magnitude 10~^ H/P in each component of the gas velocity u is generated to seed the MRI, where 
P = 2ti/VLk is the orbital period at the center of the shearing box. 



^Note that there exists a factor of \/2 difference between our definition of disk scale height H (Equation Q) and 
some other authors' in the hterature. 
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Our assumption of a net vertical magnetic field allows us to compare our models directly to 
the unstratified models presented in YMM09. The strength of this field controls the strength of the 
MRI-driven turbulence, so imposing it is an effective way to adjust the viscous stre ss generated by 



the MRI. A s imilar procedure has bee n used in other recen t publ icatioi is including I Johansen et al 



pool |200g|), iNelson k Gressell (|2010|), ISuzuki k Inutsukal (|2009l ). and ISuzuki et al.1 (|2010|). Some 
justification for this assumption can be found in the argument that at high altitude where the gas 
is tenuous , the magneto-convective motions tend to turn horizontal magnetic fields into vertical 
ones (e.g.. IProctor k Weisslll982l : iHurlburt k ToomrdllQSSl : Brandenburg et al.lll995l ). Therefore, 



we take the limit that = By = at the 2;-boundary while allowing the vertical component 
{Bz) to change freely. In terms of the vector potential A, this boundary condition translates 
to dzAx = dzAy = Az = In the horizontal directions, we use the standard sheared periodic 
boundary conditions for t he magnetic field, i .e., the field is periodic after a shift compensating for 
the Keplerian shear flow (|Hawlev et al.lll995l : see also Section [2.2p . 



We explored several different vertical boundary conditions for the gas and the fleld. We found 
that outflow bound a,ry conditions for both led to n umerical instability d uring the initial growth of 
the MRI. Although ISuzuki k Inutsukal ((20091) and ISuzuki et al.l (j2010l ) did publish models using 
outflow boundary conditions, their disk lost significant amounts of mass through the boundary 
during their run. We also found that periodic boundary conditions on both gas and field led to 
random crashes after tens of orbits of saturated turbulence for reasons that remain unclear. On 
the other hand, the combination of the vertical fleld boundary condition with periodic boundary 
conditions in the gas leads to a stable numerical result for the hundreds of orbits required for u s 
to accurately measure the evolution of particle orbital properties (see also iProctor k WeissI Il982l ) . 
The horizontal boundary conditions are again sheared periodic. In this setup, the total mass is 
conserved while saturated turbulent motion can be maintained near the vertical boundaries at a 
steady level. 

Wh ile it is interesting in it s elf to study the coronal structure a nd possibly the net flow of 
the gas (jSuzuki k Inutsukall2009l : ISuzuki et al.ll20ld : iFlock et al.ll201lh . we specifically exclude the 
modeling of the coronal regions to maintain the steadiness of the magneto-rotational turbulence 
and thus improve the statistical measurement of the particle orbital properties. The gravitational 
acceleration due to the host star across the boundary does not present any numerical difficulty. 
Even though the vertical component of the gravitational acceleration changes sign when a fluid 
element is wrapped around the vertical boundary, the vertical velocity component of the element 
retains sign such that its velocity and acceleration always have the same sign after reemerging 
from the opposite boundary, resembling material ballistically falling back down onto the disk. 
Numerically, no z-derivative of the gravitational acceleration due to the host star is ever evaluated 
(see Equation ([1])), so the discontinuity of the acceleration across the vertical boundary is not an 
issue here. Turbulence properties near the disk mid-plane, wh ere the gravitational influence of the 
turbulent gas on the particles is strongest (IQishi et al.l 120071). appear insensitive to the choic e of 



vertical boundary conditions ( Johansen et al. 20091 ; Davis et al. 2010 : Guan k Gammijboill ). 



as 
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we also demonstrate in Section [4.5[ 



For large box simulations, concerns have been raised that artificial numerical diffusion varying 
with X may occur due to the r adial dependence of the shear velocity and thus corresponding 
truncation errors on a fixed grid (I Johnson et alJl2008l ). Explicit treatment of the shear advection 
may be needed to eliminate these truncation errors. In the Pencil Code, such a technique, dub bed 
shear advection by Fourier interpolation (SAFI), has been implemented (jJohansen et al.ll2009l ). In 
our work, we implement SAFI for all of our simulations to remove this undesirable numerical effect. 



2.2. Poisson Equation with Isolated Boundary Condition in z-direction 

To find the gravitational influence of the turbulent gas on the movement of solid particles, we 
need to solve the Poisson equation for the fluctuation potential ^i: 

V^^i = 47rG/>i, (5) 

where G is the gravitational constant and pi{x,y,z) = p{x,y,z) — po{z) is the density fluctuation 
with respect to the basic state of the gas stratification po{z), which is given by Equation (j4]). Note 
that we have neglected the self-gravity of the gas in Equation ([2]) on the assumption that the 
disk is gravitationally stable (see YMM09), so the solution to Equation ([5]) does not affect the gas 
dynamics. 

Although the vertical boundary conditions for the gas density are assumed to be periodic 
in the MHD equations (see Section 12. ip , we find it inappropriate to assume the same in solving 
the Poisson equation ([5]). Consider, for example, density perturbations near the top edge of the 
computational domain. If periodic boundary conditions in the vertical direction were adopted, 
particles near the mid-plane would experience the gravity of the same density perturbations at 
almost equal distance from above and below, resulting in cancellation of the gravitational force 
exerted by these perturbations on the particles. This would amount to undesirable, potentially 
large, numerical errors for calculating the gravity of density perturbations at high altitude. 

Therefore, we instead adopt isolated or open boundary conditions in the vertical direction to 
solve Equation ([5]). With these boundary conditions, we assume that there are no fluctuations in 
density outside the vertical boundary, i.e., pi{x,y,z) = for \z\ > 2H. The boundary conditions 
in the horizontal directions remain sheared periodic. 

To implement approximate isolated vertical boundary co nditions with minimal computational 



cost, we use a variation of the iHockney &: Eastwood! (Il988l ) fast algorithm to solve the Poisson 



equation (|5l), combining it with the Fourie r interpolation technique to achieve sheared periodicity 



(jBrucker et al.l 120071 : iJohansen et al.l 120071 ). The essence of this algorithm is to double the compu- 



tational domain in the z-direction by appending regions with zero density fluctuatiorI§ and apply 



''in our case, we append pi{x, y,z) — for 2H < z < 6H so the vertical domain now covers —2H < z < 6H. 
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fast Fourier transforms to this expanded domain. Although the vertical boundary conditions re- 
main periodic with Fourier transforms, the nature of gravity means that influence from the 
periodic copies, now at least 2Lz away, of the original computational domain of size onto itself 
is significantly reduced. 

We now describe the algorithm to solve Equation ([5]) in detail. After expanding the vertical 
dimension of the computational domain and initializing the expanded region with zero density 
fluctuations, pi{x,y,z) is Fourier transformed in the y-direction into pi{x,ky, z), where ky is the 
t/-wavenumber. The result is phase shifted to recover periodicity in x-direction: pi{x,ky,z) = 
pi{x,ky, z) exp[—i{3/2)Q,Kkyx5t], where 5t is the time step. Then pi{x,ky,z) is Fourier trans- 
formed in the x-direction into pi{kx, ky, z), where k^ is the x-wavenumber. For the modes k'^ + ky > 
0, pi{kx,ky, z) is Fourier transformed in 2;-direction to find pi{kx,ky,kz), where kz is the z- 
wavenumber. The fluctuation potential in Fourier space ^i{kx, ky, kz) is calculated using the usual 
convolution theorem: 



^i(kx,ky,kz) = n Pi(kx,ky,kz), ior k"^. + k'f. > 0. (6) 



It is then inverse Fourier transformed in kz to find ^i{kx, ky, z). 

In principle, the mode kx = ky = can similarly be treated by Fourier transforms in the vertical 
direction and convolved by Equation ([6]) as above. However, this mode has an analytical solution 
that we can employ with marginal computational cost to improve the solution to Equation (l5])|l| It 
represents a horizontally infinite, vertically thin mass layer of constant density at a given altitude z, 
where the density is equal to the horizontal average of the density fluctuation at z. (The layer can 
be considered as either a slab of finite size Az and constant volume density pi{kx = 0,ky = 0,z) 
or as an infinitely thin sheet of surface density pi{kx = 0,ky = 0,z)Az; the solutions outside the 
layer are the same by symmetry and Gauss' law.) The gravitational acceleration due to this mode 
is constant above and below z. Therefore, the resultant fiuctuation potential can be calculated by 

*,fe = 0,., = 0,.)=2.G/p.fe = 0,.„=0,.')|.-.W. (7) 

where we have arbitrarily defined the reference potential to be zero at each altitude. The discretized 
version of Equation ([7]) reads 

^likx = 0,ky = 0, zj) = 27rGj2pi(.kx = 0,ky = 0, Zk)\j - k\Az'^. (8) 

k 

Equation ([8]) remains exact for this mode given that we do not consider reconstruction of the density 
field within each cell. Finally, we reverse the process to inverse Fourier transform ^i{kx, ky, z) in kx 



*We note that use of Equation ((6]) would yield a 50% relative error in the vertical acceleration due to a horizontal 
uniform slab at a distance Lz/2 away (after domain doubling) when compared to the exact acceleration derived from 
Equation 
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and ky to derive the fluctuation potential in real space y, z), incorporating the corresponding 

phase shifts with opposite sign. 

In practice, it is not necessary to allocate storage space for the whole extended domain in the 
z-direction. The appended zero density is only involved in the calculation between the forward and 
inverse Fourier transforms in x and y. In addition, the convolutions in z (Equations ([6]) and ([8])) 
for diff^erent modes {kx,ky) are independent of each other. Therefore, these convolutions can be 
distributed among processors and performed sequentially by allocating only a single one-dimensional 
working array along the z-direction. 



2.3. Particle Dynamics 

We continue to use the zero-mass approximation to model our solid particles as in YMM09. 
In this approximation, particles behaves as test particles and only respond to the gravity of the 
host star and the gas. We ignore the viscous drag forces betw een particles and gas. This remains 



a good approximation for kilometer-sized planetesimals (e.g.. iQishi et al.ll2007l ). We also ignore 
mutual gravitational interactions between particles, which helps us isolate the net effect induced 
by hydromagnetic turbulence. 

The equations of motion for each particle, therefore, become 

dxp 3 



dt "2 
~di 



^KXpY, (9a) 
^"p - ^ 2nKUp^yX - ^nKUp^xf -nj^ZpZ^ + go- V^i. (9b) 



The vector Xp = (xp, yp, Zp) is the position of the particle in the shearing box, while Up = 
{up^xjUp^y,Up^z) is the velocity of the particle relative to the background shear flow. In Equa- 
tion (I9bj) . the terms inside the parentheses stem from the linearized gravity of the host star and the 
Coriolis and centrifugal forces in the co-rotating frame of the shearing box. The remaining terms 
on the right-hand side represent the acceleration attributed to the gas: go is the gravitational 
acceleration due to the basic state of the gas stratification and has the analytical expression 

go{z) = -27T^/^GpmH erf (|:) z, (10) 

while <I?i is the gravitational potential due to density fluctuations with respect to the basic state 
and is the solution of the Poisson Equation ([5]) (Section l2.2p . Note that by separating the density 
perturbation from the equilibrium density stratiflcation and treating the gravity of the equilibrium 
state exactly, we improve the accuracy of the gravitational potential resulting from the turbulence. 

The position Xp and velocity Up of each particle is evolved in time by integrating the equations of 
motion ([9]) simultaneously with the third-order Runge-Kutta steps advancing the MHD equations. 
In addition to the Courant conditions set by the MHD equations, the time-step is limited by 
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the absolute maximum of velocity Up so that no particles can cross more than half a grid zone 
each timestep. We calculate the gradient of the fluctuation potential V<I>i on the grid and then 
quadratically interpolate it onto each particle. 

We uniformly distribute 128^ particles in a horizontal plane. We do not allow these particles 
to move until after 20 orbital periods, when the turbulence has reached a statistically steady state. 
The particles can have an initial eccentricity eo by setting an initial velocity Up so that they are at 
the apogee of their orbits (see YMM09) . They can also have an initial inclination iq by placing the 
initial particle plane at a distance to the mid-plane. Particles are wrapped around when crossing 
any of the six boundary planes of the shearing box. 



3. PROPERTIES OF MAGNETO-ROTATIONAL TURBULENCE 

In this section, we present several properties of the magneto-rotational turbulence in our 
numerical models and discuss their convergence with the box size and the resolution. 

Figure [J shows the horizontally averaged density fluctuation, inverse plasma /3, and a parame- 
ter in the mid-plane {z = 0) as a function of time, for resolutions up to 64 points per H and horizon- 
tal box dimensions up to Lx = Ly = 16 H. The magnitude of the density fluctuation is indicated by 
the rms value relative to the equilibrium density (pl)^^'^ /po, the plasma /3 = 2p.QpQcl/{\B + Sextl^) 
is th e ratio of the equili brium thermal pressure to the magnetic pressure, and a is calculated by 



(e.g., iBrandenbur g 1 19981 ) 

_ \/2 {pu^Uy - B-jcBy/po) 
3 pocf 

which includes both Reynolds and magnetic stresses^! 

The magneto-rotational turbulence in the mid-plane of our models reaches a statistically steady 
state at about t = 20P. In this saturated state, we see in general that the larger the horizontal size of 
the box, the smaller the amplitude of oscillation in the turbulence properties. Little difference exists 
between an 8x8x4H box and a 16xl6x4ff box at the same resolution, indicating convergence 
with horizontal box dimension. On the other hand, the amplitude of oscillation at saturation level 
increases with resolution for fixed box dimensions. There exists some evidence that the magnetic 
activity and thus the a value increase with resolution, especially for the case of 2x2x4// box, but 
this effect seems to become smaller with larger box sizes. We also note that there exists a curious 
trend of increase in density perturbation with time for the case of an 8x8x4H box at the resolution 
of 64 points/i/. In spite of these trends, the time average of each turbulence properties remains 
roughly the same for different resolutions and thus shows convergence. We note that the turbulence 
properties at the saturated state in the mid-plane of our vertically stratified disks agrees with those 



^The bracket ( ) denotes the horizontal average of the enclosed variable at any given altitude z and thus is a 
function of z. 
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in our unstratified disks of YMM09 (see also iGuan &: Gammid 120111 ) . The means and standard 
deviations of these properties for each set of resolution and box dimensions are reported in Table [H 
which includes the magnetic stress and the Reynolds stress as well. 

Figure [2] shows the time-averaged vertical profiles of density perturbation, inverse plasma /3, and 
OL parameter at saturation level over a period of lOOP. Each of the three properties increases with | z| , 
indicating stronger turbulence at higher altitude. The increasing activity with altitude is related 
to the increasing importance of the external magnetic field, which can be q uantified by 0e^t(z) 
(Sectionl2.1l). This a grees with the trend found in models of unstratified disks (jHawlev et al.lll995l : 
Johansen etallbood : YMM09). Notice that /3ext(±2ii') = 1.1 x 10^ > 1 at our highest altitude, and 



thus the corona is not modeled in our simulations. Given that our vertical height only amounts 
to 2H, the gas flow in our computational domain remains m arginally stable against magnetic 
buoyancy during all stag es of the development of the MRI (e.g., IShi et al.ll2010l : iGuan fc Gammie 
2011; Simon et al.ll20lll ). As shown in Figure [2l the time-averaged inverse plasma beta near the 
vertical boundary amounts to only about 0.4, and thus on average, thermal pressure still dominates 
magnetic pressure there. 



4. ORBITAL PROPERTIES OF MASSLESS PARTICLES 

We now report the response of zero-mass particles to the gravity of the density fluctuations 
of the statistically steady, numerically convergent magneto-rotational turbulence described in Sec- 
tion [3j The reference time t = in the following discussion is the time at which the turbulent gas 
reaches its saturated state and the particles are allowed to move0 



4.1. Mean Orbital Radius 



The evolution of the mean orbital radius of one particle can be found by averaging the radial 
position X of the particle over each epicycle motion. For the case of ideal unstratified disks, the 
distribution of particles in terms of the orbital radius change can be described by a time-dependent 
normal distribution centered at zero: 



f{Ax,t) 



1 



■ exp 



Ax 



2 1 



(12) 



where Ax is the orbital radius change from the initial radius at t = and ax{t) is the time- 
dependent standard deviation. We reported in YMM09 that cj^(t) depends on the properties of the 
gas disk and can be concisely expressed by 



O-x(t) 



1/2 



(13) 



®We choose this reference time to be t = 20P as reported in Section [3] See Figure [T] 
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where Cx is a dimensionless proportionality constan10and ^ = AirGpoP"^ is a dimensionless quantity 
indicating the strength of the gas gravity. Equations (fT2]) and (fT3|) demonstrate the diffusive nature 
of particle radial migration driven by magneto-rotational turbulence. 

For the case of ideal stratified disks studied here, we again find the evolution of particle orbital 
radius can be described by Equations (I12p and (I13p . The dimensionless quantity ^ is now defined 
by 

e = 47rGp„p2 ^ 4(2vr)3/2/Q^ (14) 

in terms of the mid-plane gas density pm, where Qg is the Toomre stability parameter for the 
gas. In Table [21 we report the measured values of the constant Cx from our simulations at different 
resolutions and horizontal box sizes for the case of particles with zero initial inclination. Comparing 
the value from the 2x2x4// stratified model at a resolution of 64 points/i/ with that from the 
unstratified model at the same resolution and horizontal box size (see Equation (15) of YMM09), 
we find excellent agreement in the two Cx values. 

As can be seen from Table [2] and Figure [3l the value of Cx remains roughly the same with 
different resolutions for fixed box dimensions (except the anomaly shown by the 2x2x4iJ box at 
a resolution of 64 points/ff). Conversely, it is significantly dependent on the horizontal box size 
Lfi = Lx = Ly. The larger the size, the stronger the infiuence of the turbulence on the particle 
orbital radius. This relationship can be represented by the following power-law fit: 

Cx 6.6 X 10-5 {Lh/Hf-^\ (15) 

We find no evidence of convergence with horizontal box size up to = 16H, the largest size we 
have investigated. 



Equation ()13p can be transformed into the diffusion coefficient D{J) introduced by lJohnson et al 



(|2006l ) for describing the radial random walks of orbiting particles induced by turbulent torques, 
in terms of the Keplerian orbital angular momentum J. The reader is referred to YMM09 for a 
detailed description of this procedure. We emphasize here that this transformation involves no 
assumption about the correlation time of the stochastic torques and is thus a direct measu r emen t 



of D{J). Using a heuristic choice of dimension for the diffusion coefficient, iJohnson et al.l (|2006l ) 
defined a dimensionless parameter e to represent the magnitude oi D{J). We report our measured 
values of e in Table [2j The dependence of e on the horizontal box size L/^ in our models for > AH 
is shown in Figure U] and can be written as 



e ~ 6.5 X 10-" {Lh/Hf-^^. (16) 



'^Cx as well as other dimensionless constants introduced in the following discussions may depend on the a param- 
eter. This dependency is not investigated in this work. 
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4.2. Eccentricity 

The amplitude of each epicychc oscillation of a particle gives the instantaneous eccentricity 
of the particle orbit. We reported in YMM09 that in ideal unstratified disks, the distribution of 
particles in terms of the eccentricity deviation should be a time-dependent normal distribution 
centered at zero, as long as the particles have non-negligible initial eccentricity (c.f. Equation ()12p ): 



(17) 



where Ae is the eccentricity deviation from the initial eccentricity cq at i = and the time- 
dependent standard deviation cTe(t) can be written as 



^e{t) = (^-J (^-J , (18) 

where Cg is a dimensionless proportionality constant. We emphasize that the time-dependent 
Rayleigh distributior[§ found for particles with zero (or negligible) initial eccentricity 



e2 



2al{t) 



(19) 



is a manifestation of Equation (jlTh since the eccentricity e is a positive definite quantity. Equa- 
tions ([T71) and ([19]) share the same time-dependent parameter (Te(t), and there exists no evidence 
that (Je{t) depends on the initial eccentricity cq- 

We find that the same evolution of particle distribution in eccentricity deviation also holds for 
the case of ideal stratified disks. The measured values of the constant Cg for different resolutions 
and horizontal box sizes when particles have zero initial inclination are listed in Table [2l The same 
comparison performed in Section [4.11 for the orbital radius evolution indicates that the eccentricity 
evolution in a stratified disk again agrees very well with that in an unstratified disk (see Equa- 
tion (17) of YMM09). These comparisons demonstrate that a local unstratified disk model gives 
results consistent with the mid-plane of a local stratified disk model that has the same physical 
conditions except vertical stratification. 

As in the case of orbital radius discussed in Section 14.11 the eccentricity evolution of the 
particles does not noticeably depend on the resolution for given box dimensions, while it is strongly 
affected by the horizontal box size (Table [2] and Figure [3]) . A power-law regression gives 

Ce ~ 7.2 X 10-5 {Lh/Hf-^^, (20) 

which is close to a linear relation to the horizontal box size. We discuss the box-size effect on both 
the orbital radius and the eccentricity in Section [6l 



*Note that the standard deviation of a Rayleigh distribution is equal to o"e v^(4 — n) /2. 
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The eccentricity driven by magneto-rotational turbulence enhance s orbital c r ossin g among 
planetesimals and thus increases the chance of collisions between them. Ilda et al.l (|2008l ) defined 
a dimensionless parameter 7 to represent the strength of this effect. Equation (jlSp can be used 
to estimate the value of 7, and the reader is referred to YMM09 for a detailed description of this 
procedure. We report in Table [2] our measured values of 7. As shown in Figure [H the dependence 
of 7 on the horizontal box size Lh in our models for Lh ^ 4:H can be described by 



7 ~ 3.6 X 10"^ (Lh/H) 



1.08 



(21) 



4.3. Inclination 

The only orbital property of a particle that cannot be measured in an unstratified disk model 
is the inclination i. In a stratified disk, vertical linear gravity from the host star (the third term 
in the parentheses on the right-hand side of Equation (I9b[) ) provides a restoring force such that 
particles oscillate in the z-direction about the mid-plane, as demonstrated in Figure Note 
that because the particles also experience the gravity of the gas, the period of the oscillation is 
P(l + ^/47r^) ^'^^ in the linear limit (derived from Equations ()9bp and (jlOp ). which is slightly 
shorter than the orbital period P. We can calculate the induced inclination (in radians) for a single 
particle hy i ^ (^^max — -^min) /2-R, in which z^ax and Zmin are the maximum and the minimum 
vertical positions in one oscillation, respectively, provided that z/R <^ 1. 

Figure El^a) shows the histograms of particles with zero initial inclination sorted in bins of in- 
stantaneous inclination at three different times. Similar to the eccentricity distribution for particles 
with eo = described by Equation (|19p . the inclination distribution resembles a time-dependent 
Rayleigh distribution 

/(M) = ;;^exp[-^], (22) 

where cjj(t) is a time-dependent parameter denoting the width of the distribution. 

One might expect the inclination distribution of particles with nonzero initial inclination would 
evolve as a time-dependent normal distribution with fixed center and increasing width, similar to 
the eccentricity distribution of particles with nonzero initial eccentricity. However, Figure [6]^b) 
shows different behavior: the distribution increasingly deviates from a normal distribution with 
time. The peak leans toward the mid-plane and asymmetry develops with more particles on the 
left side of the peak (i.e., smaller inclination). This indicates the diffusion is stronger near the 
mid-plane, which in fact is consistent with vertical stratification of the gas density. It would be 
interesting to measure the dependence of the diffusion coefficient on vertical height, but given 
our limited sample of initial inclinations, this remains to be investigated. We therefore focus our 
attention on the case of particles with zero initial inclination. 

Figure [7] shows ai (t) for disks with varying gravity parameter ^ and particles with varying 
initial eccentricity cq but zero initial inclination zq = 0. Note that when cri{t) is normalized by 
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S,H/R, all the curves roughly coincide. This indicates that cri{t) is linearly dependent on both 
(,H/R, but independent of cq, similar to what we reported for the evolution of the mean radius and 
the eccentricity. The results can thus be summarized by the following expression: 



where Cj is a dimensionless proportionality constant. Table [2] lists our measured values of the 
constant Ci for different resolutions and box dimensions when particles have zero initial inclination. 

We note that in contrast to the orbital radius and the eccentricity, the evolution of the orbital 
inclination is not significantly affected by the horizontal box size. The inclination of the particle 
orbits is mainly affected by the vertical structures of the density field, which in turn is related to 
the vertical scale height of the density stratification. The horizontal structures in the magneto- 
rotational turbulence may have little correlation with the vertical structures so that the horizontal 
box size has little effect on the vertical motions of the particles. 



Finally, all three components of the velocity dispersion among the particles as a function 
of time can be measured in our stratified disk models. Similarly to what we reported in YMM09, 
each component assumes the same form: 



where the index i is either x, y, or z and Si is the corresponding dimensionless proportionality 
constant. We emphasize that Sy ^ Sx/2 always holds because of the fixed ratio of amplitudes in x 
and y directions in the epicycle motions of the particles, which has been verified in our simulations. 
The values of Sx and Sz at different resolutions and box dimensions for the case of particles with 
zero initial inclination are listed in Table [2j Note that the value of Sx for the 2 x 2 x 4:H box at a 
resolution of 64 points/-ff is consistent with that from the unstratified disk with the same horizontal 
size and resolution we reported previously (see Equation (18) of YMM09), a further confirmation 
of the consistency between stratified and unstratified models discussed in Sections 14.11 and 14.21 

As shown by Figure [H we do not see evident dependence of velocity dispersion on resolution of 
our models. On the other hand, the horizontal component Sx (and thus Sy) significantly depends 
on the horizontal box size while the vertical component Sz does not. This is in accordance with 
the dependence of the three orbital properties found in previous subsections. From our measured 
values for Lh > 4i?, we quantify Sx and Sz with the following expressions: 




(23) 



4.4. Velocity Dispersion 




(24) 




(25) 
(26) 



Sz 



~ 2.8 X 10 



,-4 
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4.5. Boundary Condition Effects 

In order to confirm that our numerically convenient but perhaps physically inconsistent as- 
sumption of vertical field and periodic gas boundary conditions in the z-direction does not change 
our results significantly, we compare turbulence properties and particle dynamics from models with 
different vertical boundary conditions for the magnetic fields run at a lower resolution of 16 points 
per disk scale height H in large boxes with size 16 x 16 x 4i7. The gas density and velocity remain 
periodic. Figure [9] shows that within about one disk scale height, the vertical structure of the rms 
density perturbation, the plasma /3, and the a parameter are in quantitatively close agreement 
between the models. Note that the model with periodic boundary conditions creates a coronal 
region as low as \z\ ~ 1.8/7, where /? ~ 1. More importantly. Figure [TOl shows from top to bottom 
the time evolution in the standard deviations of particle orbital radius, eccentricity and inclination 
for the two models. The orbital radius and eccentricity agree very well irrespective of the vertical 
boundary conditions, and the effect on inclination differs by only ~30%, likely due to increased 
activity at high altitudes for the model with periodic magnetic fields. Thus, we believe that our 
results on particle dynamics are robust. 

5. PROPERTIES OF THE STOCHASTIC FORCE EXERTED BY THE GAS 

We have demonstrated in Section [3] that various properties of the saturated magneto-rotational 
turbulence in our simulations converge with both resolution and box dimensions. However, while 
converging with resolution, the response of massless particles to the gravity of the turbulent gas 
does not similarly converge with the horizontal box size, as presented in Section HI This result raises 
serious questions about the validity of using the local-shearing-box approximation to simulate the 
dynamics of any particles under gravitational influence of magneto-rotational turbulence. In order 
to understand the lack of convergence with box dimension in the orbital evolution of massless 
particles, we further study the properties of the gravitational force exerted by the turbulent gas. 

5.1. Magnitude and Fourier Amplitudes of the Force 

Table [3] lists the radial and the azimuthal components of the gravitational force exerted by the 
turbulent gas, and Fy, respectively. Since the force is stochastic, we ensemble average the root 
mean square of each force component over all particles, and then report it as time average and la 
variation scaled by (^mpHP, where nip is the mass of a particle. Note that in our models, particles 
only follow the gravitational potential generated by the gas and the host star, so the acceleration of 
each particle is independent of its mass. This is an advantage in that the measurements reported 
in this paper can be used to gauge the gravitational influence of magneto-rotational turbulence 
on solid particles or planetary objects of any mass and thus its significance compared with other 
relevant interactions. 



- 16 - 



As can be seen in Table [3l both the radial and the azimuthal components of the force generally 
increase with the horizontal box size. It is not clear if there exists any convergence in the radial 
component towards larger box size, since the 16xl6x4i? box at a resolution of 32 points per H is 
the only model that contradicts the general increasing trend. Certainly, the azimuthal component 
(i.e., the torque) shows no sign of convergence u p to the largest box siz e we have investigated, 
with a regression power-law index of about 0.6. iNelson Gressell (|201Cl ) also reported in their 
recent local models a slight trend of increasing torque with increasing horizontal box size (see their 
Figure 8 and the corresponding discussion). 

Since the density structure of the turbulent gas determines its gravitational force on the parti- 
cles, it is informative to consider how the structure changes with resolution and horizontal box size. 
The density structure of the gas is best analyzed in Fourier space, using the Fourier amplitudes 
pi {kx,ky,z) defined in Section \T2^ Figure [TT] plots the time-averaged Fourier amplitudes of the 
gas density in the mid-plane pi {kx,ky,z = 0) along either the radial or the azimuthal direction. 
The top-left panel shows the amplitude as a function of k^ for ky = 0. In general, the largest ampli- 
tude resides at the lon gest wavelength while monotonically decreasing with increasing wavenumber 



(jJohansen et al.ll2009l ). For any given box dimensions, increasing resolution has little effect on each 
amplitude, apart from extending the profile toward higher wavenumber. Flattening of the profile 
for wavelengths ^8H is hinted for 16x16x4/7 boxes. We find no evidence of a downward turn 
toward even longer wavelength. On the other hand, increasing the horizontal box size lowers 
the overall profile of the Fourier amplitudes. Interestingly, the amplitude of the longest wavelength 
mode kx = kQ = 27: /Lh remains roughly constant. 

The bottom-left panel of Figure [11] plots the summation of Fourier amplitudes over all ky at 
any given kx- In contrast to the amplitudes for ky = 0, the horizontal box size has little effect 
on the amplitudes for most of the wavelengths. Nevertheless, increasing resolution indeed increases 
the inertial range and resolves more power toward shorter wavelength. 

The right-hand column of Figure [TT] plots the azimuthal counterpart of the left-hand column. 
We find similar features in the azimuthal amplitudes as those found in the radial ones. The only 
difference is that the flat part of the spectrum in the azimuthal direction is short compared with that 
in the radial direction, as can be seen in the bottom-right panel. Power that is not captured in the 
long- wavelength range of the spectrum in the 2x2x4/7 boxes is probably the cause for artificially 
higher amplitude at {kx,ky) = (0, /cq) (the top-right panel), which may in turn be responsible for 
the anomaly found in the particle dynamics shown in Section U] 

The Fourier amplitudes shown in Figure [TT] help us understand the relationship between the 
gravitational force exerted by the gas and the horizontal box size. The radial and the azimuthal force 
components for any given (horizontal) Fourier mode k = {kx,ky) are related to the corresponding 



^Note that pi {kx, ky, z) is strictly periodic in sheared coordinates and the wavenumber fe is slightly different from 
that of a direct Fourier decomposition of pi{x,y,z)\ see Section IT2l 
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Fourier amplitude p{k) by (see Equation ([6])) 

F^{k) ~ 4vrG (I) \p{k)\ < oc L,, (27) 

F,(fc) ~ 4nG (I) |p(fc)| < 4vrG^^^ oc L,, (28) 

where we have used the facts that the dominant modes in the radial and the azimuthal directions 
are the longest wavelength ones \p{kQ,0)\ and |/5(0,A;o)|, and both are about constant irrespective 
of resolution and box size. Therefore, both components should be at most linearly proportional to 
the horizontal box size in our simulations. Although involving some simplifications, Equations ()27p 
and (I28p at least qualitatively explain the general trend of increasing stochastic force with horizontal 
box size. 



5.2. Correlation Time 

We next evaluate the correlation time Tc of the stochastic torque. This quantity is in fact not 
well defined in the literature, and different authors have different preferences for its evaluation. For 
the purpose of consistency with the frame work of radial diffusi ve migration of particles driven by 



random torque, we adopt the definition of I Johnson et al.l (j2006l ): 



T, = DiJ)/6r^, (29) 

where D{J) is the radial diffusion coefficient as a function of the Keplerian orbital angular momen- 
tum J, 6T is the fluctuating part of the torque, and the overline denotes an ensemb le average. We 



repea t the relationship between D{J) and the dimensionless parameter e here (see I Johnson et al 



20061 . and YMM09): 



To calculate the correlation time Tc with our definition, both the diffusion coefficient D{J) 
and the magnitude of the stochastic torque are needed. As emphasized in YMM09, we have all the 
information about the particle distribution against each orbital property as a function of time, and 
thus we can directly evaluate the diffusion coefficient in the context of a diffusion equation, which 
is reported in Section [4. II and Table [2j Furthermore, we also have the ensemble-averaged rms value 
of the stochastic torque, which is reported in Section [5TT] and Table [3l Equations (p9|) and ([30]) can 
then be used to calculate the correlation time Tc, and the result is listed in Table El 

As indicated in Table [3l the correlation time of the stochastic torque Tc does not appreciably 
depend on the resolution but it is significantly affected by the horizontal box size. The larger the 
horizontal box size, the longer the correlation time of the torque; this trend was also reported by 



Nelson &: Gressell (j2010l ). Up to the largest box size we have investigated (L/j < 16-ff), it is not clear 
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if Tc will converge with L/^. This behavior is similar to that of the torque magnitude, as discussed 
in Section 15. 1[ 

As a side note, one may want to evaluate the correlation time of the stochastic torque without 
referring to stochastic particle orbital evolution and invert Equation (I29p to estimate the radial 
diffusion coefficient D{J) indirectly. We have shown in YMM09, for the case of ideal unstratified 
disks, that the following formula, motivated by the definition oi D{J) under the framework of the 
Fokker-Planck formalism, serves this purpose: 



JTACFMdr 

2ACF(0) 

where ACF(r) is the autocorrelation function of the stochastic torque as a function of the time lag 
r. Here we repeat the exercise and use Equation ([3T]) to estimate Tc for the case of ideal stratified 
disks. The result is listed in Table [3] and shows remarkable consistency with the values obtained 
directly from the evolution of particle orbital distribution. Therefore, Equation (j3ip remains a 
useful estimator of the correlation time Tc- 

To further understand the trend of increasing correlation time with increasing horizontal box 
siz e, we plot the auto c orrela tion function of the stochastic torque in Figure[T2j As also demonstrated 



bv lNelson &: Gressell ([20101), the larger the box, the less oscillatory the autocorrelation function is 
and the less pronounced is the negative contribution to the radial diffusion of the particles. This 
behavior in turn is consistent with increasing values of the correlation time and the radial diffusion 
coefficient. 

In addition to the stochastic torque, a similar analysis may be conducted for the radial com- 
ponent of the gravitational force exerted by the gas. We find that the correlation time of the radial 
force also increases with the horizontal box size. However, our models indicate a much longer corre- 
lation timescale, which may amounts to several tens of orbital periods. Given that our simulations 
only lasted for about lOOP, the magnitude of the correlation time is not well calculated and much 
longer simulations may be required to cover the full spectrum of the autocorrelation function for 
the radial force. 



Johansen et al.l (|2009l ) have discussed the dominant role of the longest wavelength, purely radial 
Fourier mode k = (A;o,0) in great detail and suggested that these kind of persistent structures are 
similar to the zonal flows found in many other astrophysical systems. As noted by the authors, these 
axisymmetric structures can not generate torques that affect the orbital radius of the embedded 
particles. While the radial forces due to the zonal flows can still affect the eccentricity of the 
particles, their long correlation time indicates that they may not be responsible for the diffusive 
evolution of the particle eccentricity, as shown in YMM09 and Section[42Lj On the other hand, the 



^ It remains possible that the stochastic radial forces due to small-scale axisymmetric density perturbations drive 
the diffusive evolution of the particle eccentricities. 
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longest wavelength, purely azimuthal Fourier mode k = (0, /cq) in the gas structure^ ^1 is significant 
enough that both the radial diffusive migration and the eccentricity evolution of the particles are 
affected by the large-scale azimuthal density perturbations and thus the box dimensions. 



6. RELATION OF BOX SIZE TO PARTICLE ORBITAL EVOLUTION 



As described in Section \5\ the magnitude and the correlation time of the gravitational force 
exerted by the turbulent gas correlates with the horizontal size of our local shearing box. This 
behavior essentially explains the correlation of the stochastic evolution of the particle orbital radius 
and eccentricity with the horizontal box size reported in Section HI These results challenge the 
validity of the local shearing box to describe large-scale density structures in magneto-rotational 
turbulence. 

By measuring the two-poi nt correlation func tions of density, velocity, and magnetic field in 
magneto-rotational turbulence, iGuan et al.l ([20091) found extended density features in contrast to 
well localized magnetic structures. The authors s uggested that the density feature s are propagating 
acoustic waves excited by the turbulence (see also Heinemann Papaloizou 2009al lbl) . It is not clear 
yet what the dissipation length scale for these waves is, which is crucial to understanding whether 
a local shearing box can capture the largest-scale structures in the density fluctuations excited by 
magneto-rotational turbulence. 

On the other hand, a global disk model may require high resolu tion in order to self-consistently 
produce large-scale structures. For example. iJohansen et al.l (|2009l ) argued that the inverse cascade 
of magnetic energy from small scales to large scales might be responsible for ultimately launching 
zonal flows. To confirm this mechanism in a global context, a model that is capable of resolving at 
least the correlation lengths in the magnetic structures might be necessary. 



6.1. Bridge to a Global Disk Model? 

Given the correlation with horizontal box size reported in this work, it is not clear how a local 
shearing box can be connected to a global disk model so that both models give consistent results on 
the particle dynamics under the gravitational influence of density fluctuations in magneto-rotational 
turbulence. 

Nevertheless, we conjecture at this point that the criterion L^/R ~ 0(1) might provide a 
physical length scale for a local shearing box. At this scale, the local-shearing-box approximation 
formally breaks down since it assumes L/R <^ 1. The curvature terms become important and this 



^^Note that this mode does not represent a single shearing wave, and the wavenumber {kx, ky) is independent of 
time; see footnote [9] and Section [2.21 
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might trigger turbulent eddies to damp propagating waves. If this conjecture could be verified, 
Equations p^ . ([20|) and (j25]) would prove useful in evaluating the gravitational influence of the 
turbulent gas on particle dynamics. Substituting R for Lh in these equations, our models predict 
that 

' ~ 6.6 X 10-5 (i7/i?)-i-35 

Ce ~ 7.2 X 10-5 (H/R)-^-^^ for a ~ lO'^. (32) 
5^ ~ 1.0 X 10-4 



Nelson &: Gressell (|20ld ) presented the most recent work on stochastic planetesimal orbital 
evolution in a global disk model. Their measurement of the orbital radius and radial velocity 
dispersion can be directly compared with ours. These authors reported in their global model 
G3, whose a ~ 0.02 was the closest to ours, that crx{t) = 7.2 x 10"^ R^Jt/P and (Ju,x{t) = 
8.2 X 10"'^ Cs ^Jt/P. In this specific model, a disk aspect raticEl of H/R w 0.071 and a strength 
of gas gravity of ^ f« 2.1 at 5 AU were adopted. In comparison. Equations (fT3|) . and ([32]) 

give cj^(t) w 6.9 x 10-^ {H/^/2)^/t/P = 3.4 x lO'^ R^/tJ P and a,,^M) ~ 3 8 x IQ -^ Cs^/t/P. Our 



results are only about a factor of two lower than those of lNelson Gressell (l2010l ). Therefore, the 
criterion ~ R for a local shearing box seems to offer some consistency between global and local 
models. 

These comparisons are not meant to be accurate and conclusive, though, since the proposed 
criterion L^/R ~ 0(1) invites significant uncertainty. Further investigation into the large-scale 
structures of magneto-rotational turbulence and the discrepancy between global and local disk 
models is still needed. 



6.2. Implications for Planet Formation and Migration 

If the criterion Lh ^ R holds. Equations (I16p and (I2ip imply 
fe ~ 6.5 X 10-6 (i7/i?)-2-69 

< ^ / , ^ 1 n« for a ~ 10 (33) 

I 7 ~ 3.6 X 10-5 (i7/i?)-i-08 ^ ^ 



Consider a region at R = 5 AU in a minimum-mass solar nebula (MMSN; lHayashilll98ll ). where 
H/R « 0.08. Equations (p3]l give e ^ 7 x and 7 6 x 10"^. This value of the e parameter 
is about one order of magnitude higher than was reported in YMM09 for the case of a 2x2x2H 
unstratified box, while the value of the 7 parameter is only about a factor of three larger. Note 
that according to Equations (j33|) . both e and 7 decrease with increasing disk aspect ratio H/R, 
which often increases with increasing radial distance R. 

We suggested in YMM09 that in a typical protoplanetary disk, the radial diffusive migration of 
protoplanets induced by magneto-rotational turbulence may be unimportant compared to secular 



^See footnote [2] 
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migration. According to I Johnson et al.l (|2006l ). for the diffusive migration to be able to dominate 
over type I migration, the e parameter should be greater than or about 0.1-1 for the MMSN. 
Our new measurement remains orders of magnitude smaller than this transition value. Therefore, 
our previous conclusions on the insignificance of planetary diffusive migration may still hold even 
though e is increased by an order of mag nitude. However ^ note that recent calculations for type 

may 



I migration in non- isothermal disks (e.g., iLyra et al 



201C 



Paardekooper et al.l 120101 . 



2011 



negate the dominance of secular migration over stochastic migration driven by magneto-rotational 
turbulence. 

We also argued in YMM09 that kilometer-sized planetesimals moving in magneto-rotational 
turbulence survive mutual collisional destruc tion, except in t he inner region of a young protoplan- 
etary disk. This was based on the results of llda et al.l ([2003) for the cases of 7 = 10-3 and 10"^. 
Since our new measurement does not fall outside this range, the same conclusion still applies. 



In contrast, iNelson &: Gressell (|2010l ) argued that the velocity dispersion of kilometer-sized 
planetesimals excited by magneto-rotational turbulence might be so large that these planetesimals 
should be erosive to each other. Before a consistent scenario for the survivability of planetesimals 
could be assembled, however, several factors remain to be accounted for. First, eccentricity damp- 
ing due to tidal interaction betw een a planetesimal and its su rrounding gas acts to circularize the 
orbits of t he planetesimals (e.g.. iGoldreich &: Tremaindll980l ). This mechanism was absent in the 
models of lNelson fc Gressei (|20im and thus their measurement of the velocity dispersion of plan- 
etesimals should be considered as an upper limitEfl On th e other ha n d, usi ng analytical arguments 
to estimate the equilibrated eccentricity of planetesimals, llda et al.l (j2008l ) took three possible ec- 
centricity damping effects into account — tidal interaction, gas drag, and inelastic collisions — and 
thus their estimate of the velocity dispersion of planetesimals might be more realistic. Even though 
significant uncert ainty was i nvolv ed in the assumed strength of turbulent excitation of eccentricity 
in the models of llda et al.l (j2008l ). it could be improved by more accurate measurement of the 7 
parameter, as provided by the present work. 



2009 
2009 



Marcus et al.ll2009l') and iniproved calculatio ns on material properties (e.g.. lLeinhardt &: Stewart 



; see, e.g.. 


AsDhaue 


Leinhardt & Stewart 



Stewart h Leinhardtll2009l : IWettlaufeijl2010l ). Given all these new developments, a reconsider- 



ation of the survivability of kilometer-sized planetesimals moving in magneto-rotational turbulence 
seems warranted. 



This is not true if mutual gravitational scattering between planetesimals dominates over turbulence-driven exci- 
tation. However, this might not be the case in a typical protoplanetary environment; see YMM09. 
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7. SUMMARY 

Directly extending our previous publication (YMM09), we continue to study massless particles 
moving under the gravitational influence of density fluctuations due to saturated magneto-rotational 
turbulence in a local, isothermal, Keplerian gas disk. We include linearized vertical gravity from 
the host star and thus vertical stratification of the gas disk. For comparison, the conditions in the 
mid-plane of the vertically stratified disks are exactly the same as those in the unstratified disks of 
YMM09. 

In order to accurately measure the gravitational effect of the turbulent gas, we separate the 
gas density p into two components: the basic state po{z) for the vertical hydrostatic equilibrium 
(Equation dH) and the density deviation pi = p — po from this basic state. We use the exact 
gravitational acceleration due to the basic state (Equation ([TO]) ) and only solve the Poisson equation 
for the gravitational potential due to the density deviation (Equation ([5])). We emphasize that since 
the Poisson equation is linear in density, this approach does not assume small density fluctuations. 
Furthermore, we implement isolated boundary conditions in the vertical direction and thus any 
density fluctuation outside the vertical computational domain is neglected. 

By imposing a weak, uniform external magnetic field, we maintain a constant level of saturated 
magneto-rotational turbulence in the disk mid-plane. Several turbulence properties demonstrate 
convergence with both resolution up to 6 4 points per disk scale height H and horizontal box size up 



to 16H. The lShakura &: SunyaevI (119731 ) a parameter in the mid- plane of our models is controlled 



at the level of ~10"^. 

However, even though the properties of the turbulent gas appear numerically convergent, the 
dynamics of massless particles moving under the gravity of this turbulent gas does not converge 
with the horizontal box size L^- The larger the horizontal box size, the stronger the gravitational 
effect of the gas on the particles. Specifically, the evolution of the orbital radius, the eccentricity, 
and the horizontal velocity dispersion of the particles is roughly linearly dependent on up to 
16H. This trend was also found in our unstratified models (YMM09), and we find consistency 
between the unstratified models and the mid-plane of the stratified models. In contrast to the 
horizontal components of the particle movement, we find that the evolution of the inclination and 
the vertical velocity dispersion is not significantly affected by L^- 

The dependence of particle dynamics on the horizontal box size can be tra ced back to the 



density structure of the gas. Consistent with the large-box models studied by iJohansen et al 



(|2009l ). the longest wavelength Fourier mode dominates the density spectrum along the radial 
direction in our models. Furthermore, we find that the longest wavelength Fourier mode in the 
azimuthal direction also strongly influences the particle dynamics, leading to the diffusive evolution 
of both the orbital radius and the eccentricity of the particles. The spectral amplitudes of these 
longest wavelength modes are roughly constant against the horizontal box size. Using a simple 
single-mode analysis, we show that the linear dependence of the particle response is a natural 
outcome of these findings for the density spectrum. 
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Correlation of particle dynamics with box size poses a major difficulty for the interpretation 
of local-shearing-box models involving gravitational physics of magneto-rotational turbulence. We 
can nevertheless conjecture that Lh ~ R, where R is the distance of the box center to the host star, 
might be a natural scale of choice for a local model to approach reality. If this conjecture holds, we 
find that our previous conclusions in YMM09 on the unimportance of radial diffusive migration for 
protoplanets as well as the survivability of kilometer-sized planetesimals under collisional destruc- 
tion may still be valid. Ultimately, high-resolution global disk models and detailed comparisons 
with large-box local models might be necessary to settle this issue. 
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Table 1. Mid-plane Properties of the Magneto-Rotational Turbulence 
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Table 2. Dynamical Properties of the Massless Particles with eo = and io = 
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Table 3. Gravitational Force Exerted by the Gas and Its Correlation Time 
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Fig. 1. — Horizontally averaged properties of the modeled magneto-rotational turbulence in the 
mid-plane as a function of time. The rows of panels from top to bottom are rms density fluctuation, 
inverse plasma /3, and a-parameter, respectively. The columns are arranged with increasing reso- 
lution from left to right. Lines of different colors denote measurements from boxes with different 
dimensions. 
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Fig. 2. — Time-averaged vertical profiles of the modeled magneto-rotational turbulence at satura- 
tion stage for different resolutions and box dimensions. The panels from top to bottom are rms 
density fluctuations, inverse plasma /3, and a parameter, respectively. 
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Fig. 3. — Dimensionless proportionality constants C^, Cg, and Ci as a function of horizontal box 
size Lh = Lx = Ly. They indicate the strength of diffusion in mean orbital radius, eccentricity, 
and inclination of massless particles moving in magneto-rotational turbulence and are defined in 
Equations p^ . (fT8|) . and (p3|) . respectively. The dotted lines are power-law fits to data points with 
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Fig. 4. — Dimensionless constants e and 7 as a function of horizontal box size L^- The constant e 



i ndica tes the strength of diffusive migration driven by turbulence and was defined by I Johnson et al 
(j2006l ) , while 7 is rel ated to the streng th of eccentricity excitations (when eo = 0) due to turbulence 
and was defined by llda et al.l (|2008l ) . The dotted lines are power-law fits to data points with 
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Fig. 5. — Vertical displacement of one example massless particle moving in a stratified, turbulent 
gas disk. The simulation box has dimensions of 16xl6x4iJ, resolution of 32 points per disk scale 
height H, and strength of disk gravity C = 1- The particle has zero initial inclination. 



- 34 - 




0.002 



0.004 



0.006 
il {HI R) 



0.008 



0.01 



0.012 




-0.015 



0.005 



Fig. 6. — Histograms of particles at three different times in bins of (a) orbital inclination i when 
the particles have zero initial inclination and (b) orbital inclination deviation Ai = z — when 
the particles have an initial inclination of iq = 0.1{H/ R). The simulation box has dimensions of 
8x8x411, resolution of 32 points///, and strength of disk gravity ^ = 1. The dotted-lines in (a) 
are the best-fits using a Rayleigh distribution, Equation ([22]) . 




Fig. 7. — Parameter cTj(t) normalized by ^H/R as a function of time t, indicating the evolution 
of the width of the inclination distribution for disks with varying strength of gravity and particles 
with varying initial eccentricity. The simulation box has dimensions of 8x8x4i7 and resolution of 
32 points/ff. 
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Fig. 8. — Dimensionless proportionality constants 5^; and Sz as a function of horizontal box size L^. 
They indicate the strength of turbulent excitations for the vertical and radial velocity dispersions 
of massless particles and are defined in Equation (|24p . The dotted line is the power-law fit to Sx 
for Lh > 4.H. 
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Fig. 9. — Vertical profiles of the magneto-rotational turbulence at saturation stage for 16 x 16 x 
4:H, 16 points/If models with different vertical boundary conditions for the magnetic fields. The 
gas density and velocity for both models are periodic. The panels from top to bottom are rms 
density fluctuations, inverse plasma /?, and a parameter, respectively. These horizontally averaged 
properties are also time averaged over 20 orbital periods, a factor of five shorter than those shown 
in Figure [21 
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Fig. 10. — Evolution of the particle orbital properties for 16 x 16 x 4H, 16 points/i7 models with 
different vertical boundary conditions for the magnetic fields. The gas density and velocity for 
both models are periodic. The particles have zero initial eccentricity and inclination. The panels 
from top to bottom indicate the standard deviations of the particle distributions in orbital radius, 
eccentricity, and inclination, respectively, as a function of time. 




Fig. 11. — Time-averaged Fourier amplitudes of the density fluctuation of saturated magneto- 
rotational turbulence in the mid-plane for each resolution and box dimensions. The amplitudes 
are averaged over a period of over lOOP. The top-left panel shows \p{kx, ky = 0, z = 0)\, while the 
bottom-left panel shows the summation of \p{kx,ky,z = 0)| over ky. Similarly, the top-right panel 
shows \p{kx = 0,ky, z = 0)\, while the bottom-right panel shows the summation of \p{kx, ky,z = 0)\ 
over kx- 
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Fig. 12. — Normalized, ensemble-averaged autocorrelation function of the stochastic torque exerted 
by the turbulent gas for each resolution and box dimensions. 



